function R = compute_residual_image(image, k)
if nargin < 2
    k = 1;
end
R = image;
for i = 1:k
    R = residual_image(R);
end
end

function R = residual_image(image)
    p_rb = predict_image(image, 1, 1);
    p_lb = predict_image(image, -1, 1);
    p_rt = predict_image(image, 1, -1);
    p_lt = predict_image(image, -1, -1);
    p_mean = (p_rb + p_lb + p_rt + p_lt)/4;
    R = image(2:end-1,2:end-1) - p_mean;
end

function x = predict_image(image, xdiff, ydiff)
% xdiff, ydiff should be  +1 / -1
    [M,N] = size(image);
    cA = image(2+ydiff:M-1+ydiff,2:N-1);
    cB = image(2:M-1,2+xdiff:N-1+xdiff);
    cC = image(2+ydiff:M-1+ydiff,2+xdiff:N-1+xdiff);
    
    maxV = max(cA, cB);
    minV = min(cA, cB);
    
    x = cA+cB-cC;
    ind = (cC <= minV);
    x(ind) = maxV(ind);
    ind = (cC >= maxV);
    x(ind) = minV(ind);
end